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Abstract 

Aronnd one second after the big bang, nentrino deconpling and e^-e“ annihilation dis¬ 
tort the Fermi-Dirac spectrnm of nentrino energies. Assnming nentrinos have masses 
and can mix, we compnte the distortions using nonequilibrium thermodynamics and 
the Boltzmann equation. The flavor behavior of neutrinos is studied during and 
following the generation of the distortion. 
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1. Introduction 

Neutrinos play a significant role both in cosmology [Q and in particle physics 
1^. In cosmology, during the radiation-dominated period in the early-Universe from 
one second after the big bang to around twenty-thousand years, neutrinos are almost 
as important as photons in driving the expansion of the Universe. Furthermore, 
nucleosynthesis and the ensuing abundances of light elements are strongly influenced 
by the number and nature of neutrinos. In particle physics, it is hoped that current 
oscillation experiments coupled with solar and atmospheric neutrino observations 
will lead to insight beyond the standard model. The basic idea is that interference 
among neutrino flavors provides a sensitive experimental probe for neutrino masses 
and mixings. 

It is natural to ask whether neutrino oscillations affect the physics of the early 
Universe. Different aspects of this question have been addressed by a number of 
authors i 1 B 0 B B 0 O. 0 0. 1^. When the temperature T of the 
Universe is several MeV and higher, interactions keep neutrinos in thermal contact 
with electrons, positrons and the primodial plasma and so oscillations cannot occur. 
If thermal equilibrium were maintained for all times, then the numbers of electron, 
muon, and tan neutrinos would be equal. This would preclude flavor oscillations 
because, for instance, for each v^. that converts into a there would be a converting 
into a Vf.. 

However, around one second when T is about 1 MeV, nonequilibrium distribu¬ 
tions develop. This effect is due to the coincidence of thermal-neutrino decoupling 
and e^-e“ annihilation. As some pairs annihilate, they reheat photons and 

other electrons and positrons. The latter in turn can interact with neutrinos via the 
weak interactions, thereby slightly reheating the neutrinos. This reheating is both 
flavor and energy dependent because the neutrinos are decoupling at this time. Due 
to exchanges, electron neutrinos are heated somewhat more than muon and tan 
neutrinos. This leads to an excess of electron neutrinos over muon and tan neutri¬ 
nos. Furthermore, higher-energy neutrinos interact more strongly than lower-energy 
neutrinos so that a relative excess of higher-energy neutrinos arises. The appearance 
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of a flavor excess means that neutrino oscillations can occur. This potentially could 
affect nucleosynthesis [Q. 

The process generating the distorted distributions is sensitive to the timing of 
events in the early Universe. If the neutrinos had decoupled well before anni¬ 

hilation, then they would have maintained Fermi-Dirac distributions. Their energies 
and momenta would simply have been redshifted by the inverse scale factor R~^{t) of 
the Friedmann-Robertson-Walker (FRW) cosmology. For ultra-relativistic particles, 
such redshifts maintain standard statistical distributions. If instead the neutrinos 
had decoupled well after e’''-e“ annihilation, then they would have remained in ther¬ 
mal equilibrium during the heating process. When they decoupled later, they would 
have again maintained a redshifted Fermi-Dirac distribution. In either case, no flavor 
asymmetry would have arisen and so no oscillations would have occurred. 

The analysis of nonequilibrium neutrino distributions is not an easy undertaking, 
even in the absence of neutrino mixing. This problem has recently been treated in two 
different approaches based on the Boltzmann equation. A detailed study allowing for 
contributions from all tree-level scattering amplitudes, performing exactly many phase 
space integrations, and using numerical methods to solve the Boltzmann equation is 
presented in ref. [^. Estimates of the effects and relatively simple approximate 
analytical formulae for the distortion of distributions are provided in ref. |]K . 


When neutrinos mix, the problem of analyzing nonequilibrium neutrino distri¬ 
butions becomes signihcantly more complicated. The complications arise not only 
from the simultaneous occurrence of production and oscillations but also from indi¬ 
rect effects modifying vacuum-oscillation behavior. These effects arise from neutrino 
interactions with the background gas of electrons, positrons, and other neutrinos. 
In particular, neutrino-neutrino interactions can strongly affect oscillations because 
neutrinos in the early Universe form a dense gas and because the effects are nonlin¬ 
ear. These interactions enter as flavor off-diagonal as well as flavor diagonal terms 
in the effective hamiltonian |^. The ensuing complications can be handled using 


Hartree-Fock-like or density-matrix formalisms [18 


In previous papers [^, we have analyzed neutrino-flavor properties in the 

early Universe assuming nonequilibrium neutrino distributions based on the approx- 
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imate analytical formulae of ref. 0- These analyses therefore disregard the role 
of neutrino oscillations during the generation of the nonequilibrium distributions, 
although various tests suggest that our qualitative results are correct. 

In the current work, we obtain the equations that describe the nonequilibrium 
distortions of neutrino distributions when mixing is present, and we present numer¬ 
ical solutions of the ensuing neutrino behavior. The remainder of this introduction 
provides a guide to the structure of the paper. 

For simplicity, we assume that mixing occurs only between two neutrinos, taken 
to be Vf. and and that the third neutrino, taken as i^r, does not participate in 
the oscillations. The analysis of nonequilibrium statistics in ref. [|^ was performed 
under the assumption that muon and tan neutrinos behave identically. However, with 
neutrino mixing, the tan neutrino cannot be disregarded because it participates in 
the generation of thermal distortions via annihilations. The equations obtained 


in ref. must therefore be generalized. This is accomplished in Sect. II. 

Section III reviews our formalism for dense-neutrino oscillations in the early Uni¬ 
verse, while Sect. IV summarizes the neutrino-oscillation behavior uncovered in refs. 
TT|, O. This summary is used in the discussion of our new results in Sections VIII 


and IX. Section V combines the results of Sects. II and III to obtain equations that 
govern the production of nonequilibrium neutrino distributions in the presence of 
mixing. 

The deviations from standard thermal distributions are generated from about 0.1 
seconds to about 1.5 seconds. We call this the production phase because the flavor- 
dependent neutrino excesses are produced during this interval. By production profile, 
we mean the nonequilibrium distortions which arise at the end of the production 
phase. The term pure production is reserved for production without neutrino mixing. 
After about 1.5 seconds, neutrino decoupling is sufficiently strong that few further 
distortions in statistical distributions are generated. We call the period after about 
1.5 seconds the oscillation phase. 

One of our goals is to understand neutrino-flavor variation during the production 
phase. This requires computer simulations because the equations are nonlinear and 
relatively complicated. Issues concerning numerical methodology are discussed in 
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Sect. VI. Some improvements are made over the previous work of ref. |]^ even for 
the pure-production case. For this reason, we have repeated the analysis of the pure- 
production case in Sect. VII. 

Another goal is to obtain the production prohle in the presence of neutrino mix¬ 
ing. Production prohles for various neutrino mixing angles and mass differences are 
presented in Sect. VIII. This section also discusses the general behavior of neutrino 
oscillations and related properties. 

Section IX analyzes the oscillation phase. The main purpose is to check the results 
of refs. [0, . The qualitative flavor properties described in [^, are conhrmed. 

Small numerical differences are uncovered, however. These can be attributed to the 
use of the approximate analytical formulae for the production prohle. Finally, we 
summarize in Sect. X. 

An important improvement of the current work over refs. B0 is that the overall 
normalization of distortions is incorporated. The earlier works explicitly avoided this 
issue by considering ratios that eliminated the normalization factor. Our current 
simulations are now relatively accurate in absolute terms. At all times, the numbers 
of excess electron, muon, and tan neutrinos per cubic volume are specihed to within 
about 25%. The uncertainties are dominated by systematic effects, discussed in Sect. 
VI. 

Throughout this paper, we work in units with k = h = c = 1. The values of various 
parameters used in our current work, e.g., the baryon-to-photon ratio r], coincide with 


those of ref. . It is assumed that the chemical potentials for neutrinos are zero, so 
that the total number of neutrinos and antineutrinos is the same. 


II. Pure Production for Three Flavors 

Under the assumption that neutrinos are massless or do not mix, the distortion of 
the neutrino Fermi-Dirac distribution induced by the combination of neutrino thermal 


decoupling and annihilation has been determined in ref. [^| using the Boltz¬ 

mann equation in an expanding FRW cosmology. For this situation, the distortions 
in the muon-neutrino and tan-neutrino distributions are equal, as are the distortions 
for antineutrinos and neutrinos of the same flavor. However, in the presence of mixing 


4 


































and neutrino oscillations, these equalities do not hold, and it is necessary to gener¬ 


alize the results of ref. [^. The purpose of this section is to obtain the generalized 
equations. 

We begin by listing four of the more important approximations made in the analy¬ 
sis. First, at the energies and temperatures of interest, the difference between Fermi- 
Dirac and Maxwell-Boltzmann statistics is unimportant because there is little fermion 
degeneracy. For calculational purposes it is convenient to use Maxwell-Boltzmann 
statistics. A second approximation is the restriction to dominant weak-interaction 
effects, of order G\. This is an excellent approximation during the time interval of 
interest in the early Universe. 

A third approximation concerns the value of 


m = T-,iT-i , 


( 2 . 1 ) 


where is the photon temperature and T is the neutrino temperature. The quantity 
6{t) is a measure of the photon-neutrino temperature difference Ty — T, to which the 
distortions of the neutrino distributions are proportional. An exact computation of 
5it) is infeasible. However, it may be approximated by where 6o(t) is obtained 
assuming that neutrinos are always thermally decoupled and that they do not share 
in the heat released by annihilations. While this approximation is poor for 

early times when the temperature is above ~ 5 MeV, the corresponding values of 
6(t) are relatively small, and so the error introduced in the distortion is minimal. A 
correction to 6o(t) partially compensating for the approximation is given in ref. [|^ . 


A fourth approximation is to set the electron mass mg to zero. This is particularly 
convenient when computing scattering amplitudes. The approximation is reasonable 
for T > 0.5 MeV. It does not introduce a large error for T < 0.5 MeV because most 
of the production occurs at temperatures above 0.5 MeV. 

We next introduce the neutrino distortions. When neutrinos are in thermal equi¬ 
librium, their energy distribution is governed by the Maxwell-Boltzmann factor fo{E) 
given by 


fo{E) = exp - 


E 


( 2 . 2 ) 
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where E is the neutrino energy. As T drops below ~ 5 MeV, neutrinos decouple 
thermally. The distribution {E,t) for each neutrino Uf of flavor /, Uf = Uf,, I'r, 
9f,, 9^ or 9^-, then deviates from the Maxwell-Boltzmann one. We write 


/„, (£,()=/„(B) + A„, (B,() 


(2.3) 


This equation dehnes the distortions {E,t). 

For early times t when the system is hot, t) = 0, fuf{E) = fo{E) and there 

is no deviation from standard statistics. For later times, {E,t) is nonzero but 
small compared to fo{E) because the neutrinos are weakly interacting and are only 
mildly sensitive to the reheating from e+-e“ annihilation. Thus, {E,t) can be 
treated as a perturbation. Since electron neutrinos are somewhat more sensitive to 
annihilation than muon or tan neutrinos, Aj,^ is larger than Aj,^ or A^^. 

The procedure for obtaining equations for the A^^^ is given in ref. 0. We therefore 
restrict ourselves here to outlining our derivation of the generalized equations for the 
distortions. Note that, throughout the calculations, each neutrino momentum p and 
energy E can be equated. This is an excellent approximation since the neutrinos are 
ultra-relativistic at the temperatures of interest. Also, in working with the Boltzmann 
equation for f^fiE) in an expanding Universe, it is convenient to use the variable E/T 
whenever possible since, for relativistic particles, it does not redshift as the Universe 
expands. The FRW scale factor R{t) evolves to maintain T{t)R{t) constant to a good 
approximation because neutrinos are decoupling and do not participate heavily in the 
reheating. In contrast, the photon temperature Ty rises relative to T during 
annihilations. 

The squared matrix elements of neutrino-scattering processes appear in the Boltz¬ 
mann equation. There are 42 different basic processes: 


I'a + e~ ^ + e~ , 

I'a + e^ ^ Ua + e^ , 

I'a + I'b ^ I'a + I'b , 


+ l^b + ^b , 

9a + e~ ^ 9a + e~ 

9a + ^ 9a + e+ 

9a + 9b^ 9a + 9i, , 
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+ , a i^b 


l^a + ^b^ 


(2.4) 


where a and b stand for e, /r, or r. In the Boltzmann eqnation, an initial nentrino 
is distingnished. This leads to 24 additional cases. Thns, there are a total of 66 
distingnished-nentrino processes. For each particnlar nentrino or antinentrino, there 
are 11 cases. 

In the calcnlation, hard scattering of nentrinos off nncleons can be neglected be- 
canse the nncleon density is small. Also, nncleons do not affect nentrino oscillations 
throngh forward scattering becanse electron and mnon nentrinos are scattered in the 
same way. As noted above, processes involving electrons and positrons tend to reheat 
the nentrinos as e"*" and e~ annihilate, and electron nentrinos participate more in the 
reheating dne to charge-cnrrent interactions from exchange. In contrast, pnre 
nentrino and antinentrino processes tend to eqnilibriate the distribntions. 

The electron and positron Maxwell-Boltzmann distribntions are 

fe± {Ee,t) = exp 

Since the electron mass is set to zero, = Pe- Expressing Eq. (|2.5|) in terms of the 
nentrino temperatnre T and the fractional temperatnre difference 6{t) gives 

/e± (Pe,t) = /o(Pe) (^1 + + • • •) • (2.6) 


The Boltzmann eqnation itself and the sqnared matrix elements for all the pro¬ 
cesses in Eq. (|2.4|) are given in ref. . To obtain eqnations for the distortions, Eqs. 
and (|2.6|) and the sqnared matrix elements are snbstitnted into the Boltzmann 


eqnation, and an expansion is performed in powers of Ap^ and 6(t). The term 
linear in these qnantities provides differential eqnations for the Aj,^. For calcnlational 
pnrposes, it is nsefnl to consider the distortions A,^^ and other qnantities as fnnctions 
of E/T and t, rather than E and t. 

The resnlting differential eqnations can be simplified by performing many of the 
phase space integrations. After some calcnlation, we obtain for the six distortions 
Ayj. a set of six linear conpled differential eqnations, given by 


dA 


Uf 


dt 


{§,t) 
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where the A, B, and C coefficients are determined. They are presented in Appendix 

A. 


Since > 0, the A terms indnce damping. The B terms indnce reheating 
of nentrinos, arising from interactions with positrons and electrons. Note that onr 
conventions for the A and B coefficients differ from ref. by a factor of E/ (Tvr^), 
inclnded here as part of the overall normalization constant in Eq. The C terms 

represent the conpling of nentrinos of different energies. This is the origin of the 
integration over the hnal state nentrino energy E'. Note that all the coefficients are 
independent of t, since they are fnnctions only of E/T and E'/T. 

We remark that Eqs. (|2.7|) correctly rednce to the corresponding eqnationsj^ of ref. 
r5| in the limits 




A. 




A,.. 


A.. 


A. 


( 2 . 8 ) 


III. Pure Oscillations for Two Flavors 

In this section, we snmmarize the formalism of onr earlier work for treating nen¬ 
trino oscillations in the early Universe. We consider oscillations between two flavors 
Pe and 

The most convenient formnlation nses the three-component vectors 
v{E, t) = (ylue - 2Jm(z/]z/^)) , 

w{E, t) = - hjp/,, 2i?e(P^P^), 2Jm(P^P^)) , (3.1) 

where Uf and Pj are helds for nentrinos and antinentrinos of flavor /, respectively. We 
nse the normalization convention that ^'jz/j(E) is the nnmber of nentrinos of flavor / 
with energy E in a comoving volnme of volnme per nnit E/T. The scaling factor 
a can be chosen arbitrarily. It increases with the redshift scale R{t) as the Universe 
expands. 

^Note that the factor (c-|-8) in Eq. (2.lid) of ref. ||T^ should be {c + 7), and the factor (6-|-c-|-14) 
in Eq. (2.14b) should be {b + c + 13). 











The oscillations equations resemble the motion of a particle in a magnetic held: 

= v{E) X B,{E) , = , (3.2) 

where the effective magnetic helds are 

By{E,t) = — Vyy — {ycp+{E) + Vcp-)ei , 

Bn,{E, t) = ^ + C - {VcP^{E) - Vcp-)e, . (3.3) 

The next few paragraphs explain the terms in Eq. (|3.3|) . 

Vacuum oscillations are produced by the term A/(2E), where 

A = A (cos 20, — sin 20, 0) . (3.4) 

Here, 0 is the vacuum mixing angle and A is the mass-squared difference of vacuum- 
mass-eigenstate neutrinos: A = m\ — m\. For pure vacuum oscillations, neutrinos 
are in mass eigenstates when v and A are aligned. 

The interaction of neutrinos with background electrons and positrons is governed 
by two potentials, Vcp+ and Vcp-, where Vcp+ is CP-conserving and Vcp- is CP- 
violating: 

Vcp+{E^ t) = —2y/2GpE (Pe- + Pe- + Pe+ + Pe+) /; 

Vcp-{t) = V^Gp {rie-— ne+) , (3.5) 

where nf,pf and pf are respectively the number density, the energy density and the 
pressure of the charged lepton /. As the Universe expands, these quantities decrease. 
The effective magnetic helds associated with Vcp+ and Vcp- are aligned along the 
1-axis hxed as ii = (1, 0, 0). When neutrinos are in pure havor eigenstates, the v are 
aligned with ei. 

Finally, provides the neutrino-neutrino interactions. It is given by 

V.(() = ((S) - («!•» , (3.6) 

where 

oo oo 

= J—v{E,t) , {w{t)) = J—w{E,t) . (3.7) 

0 0 
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The asterisk in Eqs. (^) and (|3.6|) indicates that the sign of the third component 
of a vector is reversed, e.g., { 101 , 102 , 103 )* = {wi,W2, —W3). Since renders the 
oscillation eqnations in Eq. ( |3.2|) nonlinear, we often refer to it as the nonlinear term. 

IV. Summary of Results for Pure Neutrino Oscillations 

For pnrposes of comparison with the new resnlts given in the present paper, this 
section provides a snmmary of the main featnres of nentrino behavior in the early 


Universe described in refs. [|lTl |^. Numerical solution of Eq. (|3.2|) for A > 0 and 0 < 
6* < 7r/4 reveals that neutrino oscillations exhibit a number of effects, all attributable 
to the nonlinear neutrino-neutrino interaction term. 

One feature of the neutrino behavior is its remarkable smoothness. As the Uni¬ 
verse expands, the four interaction terms in Eq. ( p.3|) change. In the absence of the 
nonlinear term, the changes induce temporary oscillatory behavior. However, when 
the nonlinear term is present, no oscillations are observed. 

Another feature of the behavior concerns decoherence. Without the nonlinear 
term, decoherence is a collective effect due to the different oscillation times of indi¬ 
vidual neutrinos. Individual neutrinos nonetheless undergo oscillatory behavior. In 
contrast, when the nonlinear term is present, individual neutrino vectors are approx¬ 
imately aligned and so do not exhibit independent oscillatory behavior. 

A third feature of the nonlinear system is that individual neutrinos maintain them¬ 
selves in configurations that are approximate nonlinear mass eigenstates (ANME). An 
instantaneous nonlinear mass eigenstate diagonalizes the hamiltonian at a given time. 
This occurs for an individual neutrino when its vector and magnetic field point in 


the same direction, so the right-hand side of Eq. (^) vanishes and the change in 
flavor content is momentarily zero. As the Universe expands, magnetic fields evolve. 
The ANME property means that a neutrino vector tends to follow changes in the 
associated magnetic field. It can be shown that alignment is a consequence of the 
ANME property and the dominance of the neutrino-neutrino interaction. 

A fourth feature of neutrino behavior is CP suppression. Although CP-violating 
interactions are present, the buildup of asymmetry between neutrinos and antineu¬ 
trinos is several magnitudes smaller than in the absence of the nonlinear neutrino- 
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neutrino interactions. The CP-suppression mechanism can also be understood as a 
consequence of the ANME property. It involves a partial cancellation between the 
CP-asymmetric electron-neutrino and neutrino-neutrino interactions. This cancella¬ 
tion can be used to develop an analytical approximation scheme for the average and 
individual neutrino vectors. The scheme reproduces individual neutrino vectors to a 
few percent and the average vector to better than one percent. 

The neutrino behavior also exhibits planarity, i.e., the third components of neu¬ 
trino vectors are much smaller than the first and second components. In a linear 
system for which decoherence has set in, planarity holds for the average vector but 
not for individual neutrino vectors. Since instantaneous nonlinear mass eigenstates 
have zero third components, the ANME property implies that the third components 
of individual vectors are also small for the full nonlinear system. 

Since the third components are small, it is tempting to conclude that they are 
unimportant. This is false. Disregarding them, which is equivalent to disregarding 
the imaginary part of the off-diagonal flavor term, is a poor approximation. While 
the first two components cancel to a large extent in the nonlinear term of Eq. ( |3.6|) , 
the third components tend to add because of the asterisk operation. Numerically, we 
find the third component of Wi/ is of the same order of magnitude as the first two 
components. 

Another feature of the neutrino behavior can be traced directly to the third com¬ 
ponent of Vvy. When A < 10“® eV^, ANME configurations initially point approxi¬ 
mately along the 1-axis. As the Universe expands and energies redshift, the vacuum 
term grows, becoming dominant some tens of seconds later. The magnetic fields 
then point toward A, so ANME conhgurations are close to vacuum-mass eigenstates. 
Surprisingly, however, the neutrino vectors rotate toward vacuum-mass eigenstates 
earlier than expected. This effect, called precocious rotation, arises because a nega¬ 
tive contribution develops to the third component of Wi/, which causes a rotation of 
the neutrino vectors in the 1-2 plane. 

In passing, we mention two other effects involving the imaginary part of the off- 


diagonal flavor term, arising in the A < 0 parameter region |^. First, there exists 
a flavor-conversion mechanism that is completely different from the MSW effect |P0||. 
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It is efficient and occurs even for very small mixing angles. The second effect, called 
self-maintained coherence, is a collective mode of the nonlinear system in which a 
large fraction of neutrinos oscillate in unison. The effect in the pure neutrino gas 
0 differs from that in the mixed neutrino-antineutrino gas [|T^ , where the third 


component of Vyy plays an important role. 


V. Equations for Combined Production and Oscillations 

If neutrino masses and mixings are nonzero, neutrino oscillations begin as dis¬ 
tortions develop in the neutrino distributions. In principle, the oscillations could 
signihcantly affect the flavor content when the photon reheating and neutrino decou¬ 
pling are complete. For large mixing angle, for instance, any excess electron neutrinos 
produced could oscillate into muon neutrinos, potentially producing a muon-neutrino 
excess that in the absence of oscillations would not arise. Furthermore, any change 
in flavor content due to oscillations affects the rates of the scattering processes in the 
Boltzmann equation. The goal of this section is to obtain equations governing the 
production of nonequilibrium distributions when neutrino mixing is present. 

During an inhnitestimal time interval, the change in the neutrino distributions is 
due to the change created by the Boltzmann equation plus the change due to neutrino 
oscillations. Hence, it suffices to add the contributions from pure production and from 
pure oscillations. This involves combining the equations in Sects. II and III. However, 
the formulae in Sect. H are expressed in a different basis from those in Sect. HI. For 
numerical purposes, we elect to use the vector-type variables of Sect. HI. We therefore 
begin this section by rewriting the pure-production equations of Sect. H. 

The hrst step in this process is to relate the distortions to the density-like 
variables The connection is given by 

. (5.1) 

Note that the factor E'^a^T / (27r^) is independent of time. Multiplying both sides of 
Eq. (E3) by this factor leads to the reformulation of the pure-production equations 
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as 


dull 


AGlT^ 

71 ' 




+ B„,(f)i(*) + E/^C'w(f'f)4''/'(^''‘)) • (8-2) 


''f' 0 


where the tilde coefficients are defined by 


B _ p (®^) f^\ ^ (E E'\ — fi (E E'\ /c o', 

J^Uf yrpj — JDyj 2^2 yj' J ’ \T^ T J — \T ^ T ) ' 


The A coefficients remain unchanged. Since is time independent, so is \§)- 
The second step is to convert to the vector formulation. To accompany the 1, 
2 and 3 neutrino-vector components specified in Sect. Ill, we define new 0 and 4 
components by 


^0 (f, + vluf}j {E, t) , V4 (f, t) = ulur {E, t) . (5.4) 

and 

Wo (f, t) = {ulue + {E, t) , Wl (f, t) = ulur {E, t) . (5.5) 

Here and henceforth, we use a bar over a neutrino-vector index to help distinguish 
neutrinos from antineutrinos. Also, to render the variable dependence the same for 
production and oscillations, we use the arguments E/T and t instead of E and t for 
the new components vq, V 4 , Wq and w^. The zero component vq is associated with the 
electron-muon neutrino excess in the canonical comoving volume, while the fourth 
component V 4 is the tan-neutrino excess in this volume. Corresponding statements 
hold for the antineutrinos. 


When reexpressed in vector notation, Eq. ([^.7|) for pure production becomes 



where a = 0,1, or 3, 5 = 0,1, or 3 and 6 = 0,1, or 3. The new coefficients Aab, 
Ba, Cab, and Cai in this equation are linear combinations of the previously defined 
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coefficients and They are given in Appendix B. The new coefficients 

are fnnctions of nentrino energy and temperatnre, but again only in the combinations 
E/T and E'/T. They appear due to the change of basis from the vectors Ci and 62 
to the vectors ei + 62 and 61 — 62 . 

A similar derivation holds for antineutrinos, yielding 


^ (f > E ^-ab (f) ^b (f > t) + B-a (f) 5 {t) 

00 j 

60 b 

The coefficients and Bg, Cab and Cgi are also provided in Appendix B. Note that 
if the hrst index of a coefficient is unbarred then it enters in a neutrino equation. 
If instead the hrst index is barred, then the coefficient enters in an antineutrino 
equation. Equations and (|5.7|) complete the change of basis to vector notation 
for the pure-pro duct ion case discussed in Sect. II. 

The equations for the combined system with simultaneous production and oscil¬ 
lations follow from the above analysis. For neutrinos, adding the right-hand sides of 
Eqs. & and ( |5.6|) gives 



dvg 

dt 


(At) = (f;(f,t) 



a 


7J-3 


E ^'^b (f) Vb (f, t) -h Bg (f) S (t) 


y 7 (f. f) V, (f. t) + Y.J (f. f) (f. () 


b 0 


For antineutrinos, adding the right-hand sides of Eqs. (|3.2|) and gives 
dwa 


dt 


(f ,f) = (w (f ,f) X B^y^ 


AGlT^ 


TT"’ 


E ^ab (f ) (f, t) + Bg (f) 5 (f) 


E (f. f) (f. *) + E 7 (I, f) (S, () 


b 0 


(6.8) 


(5.9) 
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Equations (|5.8|) and (|5.9|) are the desired equations. Note that in Eqs. 
we define 


and (1^ ) 


(il (E) X = 0 , for a = 0 and 4 , 

(w (E) X Bw^_ = 0 , for a = 0 and 4 . (5.10) 

The 1, 2 and 3 components are computed normally. 


VI. Technical Issues and Numerical Methodology 

In this section, we first describe the discretization of Eqs. ( b.8|) and (|5.9|) - We 
next present some analytical results for the resulting equations and discuss the im¬ 
plementation of the approximation for 6{t), already mentioned in Sect. II, along with 
a method for partially compensating the approximation of zero electron and positron 
mass. Finally, we present some details about the algorithm used and the range of 
data obtained. 

Numerical solution of Eqs. ( |5.8| ) and ( |5.9|) requires conversion of the continuous 
energy variable to a discrete one. We allow E/T to vary from (E/T)^.^ to (E/T)^^^ 
and divide this interval into Ne bins of equal spacing. Indices j and k ranging 
from 1 to Ne are used to label quantities associated with neutrino and antineutrino 
energy bins, respectively. Most simulations were performed with (E/T)^.^ = 0.1, 
iE/T)max = 20.0, and Ne = 61. 

We define discrete density-like variables associated with the jth and kth bins by 

ufuj, = [e^) , z/f z/J, = v]Vf, {e^) . (6.1) 

With this definition, represents the number of neutrinos of flavor / with energies 
between E^ and E^ -|- NE^ in the canonical comoving volume a^(t), while is 

the number of antineutrinos of flavor / in the fcth energy bin in the volume a^{t). 
Discrete neutrino and antineutrino vectors associated with the jth and fcth bins can 
similarly be introduced, via 

A JT'j /\ 

= , W^{t) = ■ ( 6 . 2 ) 

Differential equations for v{ can be obtained by multiplication of Eq. (|5.8|) with 
AE^/T and replacement of the continuous integral over E'/T by a discrete sum. 
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Since AE^/T is time independent, it can be moved inside the time derivative on the 
left-hand side. For antineutrinos a similar procedure is performed to give equations 
for w^. The resulting discrete equations resemble the continuum equations (|5.8|) and 


(5.9) with the identihcations 


E ^ E\ E' ^ E^ 





Va t) Vi{t) , Wa (^, t) ^ W^{t) , 

Ai (f) - = A,, (f) , B, (I) ^Bi= ^B, (f) , 

'^cd \ ti T ) ^Cd — \ t i T ) ’ 

B, {E, t) ^ Bi{t) = B, [e\ t) , B^ {E, t) ^ Blit) = B^ {e\ t) , 
Vcp+ (E) ^ Vl;p+ = VcpA [e^ in B, and B^ of Eq. (jp) , 


{va{t)) in of Eq. (3.6) , (6.3) 

j k 

where the label c can be a or d and d can be b or b, and where the indices i and m 
stand for j ot k. 

In the absence of production, i.e., setting = B^ = C*™ = 0, the equations 
reduce to those in ref. ||TI|, |T^ for pure oscillations. In the absence of oscillations, i.e.. 
By = By, = 0, the discrete set of differential equations can be analytically integrated, 
as we show next. 

The hrst step is to convert the time integration variable t to the neutrino temper¬ 
ature T, via 


T 


1 



(6.4) 


where 


Kq = 




(6.5) 


Here, is Newton’s constant and counts the total number of effectively massless 
degrees of freedom. For T > 1 MeV, = 10.75. 
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The structure of the discrete version of Eqs. 
becomes 


|) and or of Eq. (|^) then 


dxr (T) 
dT 


= -koT^ {BrHT)+CrsXs (T)) 


( 6 . 6 ) 


where B is related to the B coefficients and C contains both C and A coefficients. Note 
that both the sets of coefficients B and C are independent of time and temperature. 
In Eq. (|6.6| ), 


^0 =- n 

KOTT'^ 


( 6 . 7 ) 


and r, s represent both flavor and energy-binning indices: r ^ cA with c ranging 
over six values and i ranging from 1 to Ne- In the formulation of Eqs. (^.8|) and (|5.9|) , 
X represents either v or w, and c ranges over 0,1,4, 0,1,4. For the formulation of Eq. 
(E3). X represents one of the six distortions and c ranges 
The solution to Eq. 


over I/e, Z/;,, Z^r, Z^e, I'fi, 


IS 


r df T^6 (t) 

exp (^(f3-T3) c] 

JT ^ 

[ V3 ^ MJ 




( 6 , 8 ) 


Although the QNe x QNe matrix C is not symmetric, it does have a complete set 
of right eigenvectors defined by 




( 6 . 9 ) 


where q = 1,..., GNe- The are the right eigenvalues of C. The are not 
orthogonal but they are complete. Thus, the vector B can be expanded in terms of 
the using 

, ( 6 . 10 ) 

<? 

where the /9(g) are expansion coefficients. The solution (|0|) becomes 

x, (T) = fco E dTfH (f) exp (f 3 - T^) . (6.11) 

Equation (|6.11|) provides an efficient way to obtain results for the pure-production 
case. The integrals in (|6.11|) can be computed rapidly by various methods. One 
needs ReA*-'^^ < 0 for convergence. For different values of Ne-, we have computed the 
right eigenvalues of C and verified convergence for sufficiently large Ne- This analysis 
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aided our determination of the minimum size of Ne necessary for numerically accurate 
results. As a side remark, we note that a reasonable approximation to pure production 
can be obtained by saturating the sum in Eq. ( p.ll| ) with terms corresponding to those 
eigenvalues with the largest real parts. 

To integrate the distortion equations, a formula for 6 {t) is needed. We follow ref. 
r5| and approximate 6 (t) by So(t), given by 

\ 1/3 


So{t) = 


12 


- 1 


( 6 . 12 ) 


A + z^Ki{z)+Az^K2{z)^ 

where z = rrie/T^ and Ki{z) and K 2 {z) are modihed Bessel functions. Recall that 
6 o{t) is obtained under the assumption that neutrinos are thermally decoupled for all 
times. 

The direct use of Eq. (|6.12|) for arbitrarily early times leads to an inconsistency. 
It implies an early-time behavior 


1 frUe 
Oo(t) ^ — — 

ov ; 36 V T 


( 6 . 13 ) 


for T S> rUe. If this expression were used in the pure-production equations, one would 
End that the early-time behavior of 




^Vf 




(6.14) 


where ^ time-independent constant. Since the total excess density is ob¬ 

tained by multiplying by E^dE/2n‘^ and integrating over E, it would follow that the 
excess density of electron neutrinos over muon neutrinos increases at earlier times. 
If this result were correct, it would pose a severe numerical difficulty. When both 
production and oscillations are treated, the neutrino-neutrino potential would be im¬ 
portant for T < 100 MeV. This would necessitate starting the integration of Eqs. 


and (|5.9|) at T ~ 100 MeV. To integrate over the large temperature range from 
~ 100 MeV down to ~ 1 MeV would be prohibitive in computer time. 

Fortunately, the early-time behavior in Eq. (|6.13|) is incorrect. The source of the 
difficulty is the assumption that neutrinos are thermally decoupled for all times. It 
must be true that for T > 10 MeV neutrinos are in thermal equilibrium with photons. 
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so that <5 is exponentially small. The results of refs. and suggest that neutrinos 
decouple in the temperature range from 2 MeV to 5 MeV. Higher-energy neutrinos 
decouple later, at a lower temperature. Also, electron neutrinos decouple somewhat 
later than muon and tan neutrinos. In any case, it is evident that S should be set to 
zero at sufficiently early times. 

An exact formula for 5(t) is not available. Instead, we make the simple approxi¬ 
mation 

5(t) = hoW , fort >4 , 

5(t) = 0 , for t < tg , (6.15) 

where tg is a cutoff time. Numerically, t* is the time at which we start the integration 
of Eqs. (^]8D and 

In principle, a better approximation can be made by taking tg to be a function of 
bin energy. In other words, 5 is set to zero earlier for low-energy neutrinos and later 
for high-energy neutrinos. If one is interested only in results for t > 0.3 seconds, then 
the detailed manner in which 6 is turned on is unimportant. For simplicity, we take 
tg to be energy-bin independent. We denote by Tg the neutrino temperature at the 
time tg. 

The approximation (|6.15| ) is a significant improvement over the use of at all 
times. It is also better than incorporating the 6 T^/T correction to 6 provided in ref. 
T5| . For numerical purposes, we take Tg = 3.0 MeV, corresponding to a starting 


time tg of about 0.082 seconds. We have varied Tg somewhat and checked that hnal 
production results are relatively insensitive to this choice. Note that, although it is 
important to turn off 6 for early times in dealing with neutrino oscillations, it is less 
important in the pure-production case if one is interested in results at late times. 
Hence, using the approximation 5o at all times is unlikely to affect the conclusions of 
ref. concerning the primordial “^He abundance. 

Another approximation made is to set the electron and positron mass to zero. 
When the temperature drops below 0.5 MeV, electrons and positrons quickly annihi¬ 
late and production rapidly ceases. This is not correctly incorporated in the formulae 
for the scattering processes when rrie = 0 is used, so some error is made below 1 
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MeV. A simple way to compensate for this is to stop the production prematurely. 
We chose to terminate our numerical integrations involving production processes at 
the temperature Tf = 0.75 MeV, which corresponds to a time tf of about 1.312 sec¬ 


onds. Beyond Tf, we proceed with integration of the pure oscillation equations (|3.2| ) 
instead. Varying Tf from 0.9 to 0.5 MeV changes the distortion functions by about 
20%. Hence, the 20% uncertainty in selecting Ty is the biggest uncertainty in our 
results. 

To minimize the chances of programming errors, at all stages two independent 
programs were written and results were checked to double-precision machine accuracy. 
Desktop Hewlett-Packard and Sun work stations were used. 

The numerical integrations were performed with a fourth-order Runge-Kutta al¬ 
gorithm applied to Eqs. (|5.8| ) and (^.9|) . The scale factor a(t) at = 3.0 MeV was 
taken to be 0.5 MeV“^ to facilitate comparison with the results of ref. [^. The point 
is that when T* = 1.5 MeV, a{t) becomes 1.0 MeV“^, which is the value selected 


in ref. |^. Since the A, B, and C coefficients are time-independent, they could be 
computed once initially and stored in memory. 

The integration combining production and oscillation proceeds about one-hundred 
times slower than pure oscillation integration due to the increased amount of arith¬ 
metic. For this reason, verihcation of all the runs in ref. was impractical. However, 


sufficiently many simulations were performed to con&m the conclusions of ref. [|l3 
In addition to the situation without oscillations, we studied the cases sin^ 29 = 0.81 
with A = 10“^^ eV^, A = 10“® eV^, A = 10“’^ eV^, A = 10“® eV^, and A = 10“^ 
eV^, and the cases sin^ 29 = 0.25 with A = 10“^^ eV^, A = 10“® eV^, and A = 10“® 
eV2. 


VII. Results for Pure Production without Oscillations 

For the case without flavor mixing, the improvements considered in Sect. VI lead 
to somewhat different nonequilibrium distortion prohles than those of refs. |^, K . 
This section presents our results and discusses these differences. 

Figure 1 displays our hnal production prohles at Tf = 0.75 MeV for zero mixing 
angle. The curves plotted are h'lh'e{E)/T, /T, their sum, and their difference. 
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as functions of neutrino energy E. The continuous distortion densities ujufi^E) are 
defined in Eq. (|5.1| ). We remind the reader that i'ji'f{E)AE/T represents the excess 
number of neutrinos of flavor / between E and E + AE in the canonical comoving 
volume. Each curve therefore represents a neutrino excess per unit energy over the 
Maxwell-Boltzmann distribution in the volume a^, and the area under a curve is the 
total excess in that volume. 

The production prohle for the tau neutrinos is identical to that for the muon 
neutrino because, in the absence of mixing, the system is symmetric under the inter¬ 
change of muon and tau flavors. Also, the distortions for antineutrinos are the same 
as for the corresponding neutrinos because the system is CP symmetric when 6* = 0. 

The electron- and muon-neutrino curves are negative for small E, indicating a 
dehcit of neutrinos in the low-energy region compared to the Maxwell-Boltzmann 
distribution. Physically, the effect arises from the reheating of neutrinos by 
annihilation, which shifts some low-energy neutrinos to higher energies and thereby 
reduces the number in the low-energy region. As can be seen from Figure 1, the dehcit 
largely cancels in the neutrino-difference prohle, which represents the diherence be¬ 
tween electron- and muon-neutrino number densities. The havor-density nonequilib¬ 
rium distortions peak at E ^ AT, while the diherence curve peaks at E ^ 3.3T. All 
curves have exponential tails that become small for E > 13T. 

The distortion curves in Figure 1 evolve with time for t > tf due to the expansion 
of the Universe and the associated energy redshifts. For T <0.75 MeV, the curves 
maintain the same overall shape but are compressed horizontally and expanded ver¬ 
tically by the factor T//T, where Tf = 0.75 MeV. Thus, at later times the peaks at 
~ 3.0 MeV approach the origin and increase in height. The areas under the curves 
remain constant, since neutrinos are neither created nor destroyed when t > tf. 

It is useful to have analytical expressions representing the pure-production prohles 
in Figure 1. We hnd the electron- and muon-neutrino prohle curves are accurately 
reproduced by the expressions 

ulve{E)/T ^ (-1.3567 x 10“^ + 1.9692 x lO'^a; + 2.4600 x 10"V 

\ 

+5.9818 X 10"V - 7.5683 x 10“V) — exp(-a;) , 

/ 27r^ 
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(7.1) 


vlv^{E)/T ^ (-1.0568 X lO"^ + 1.0668 x 10“^x + 1.4263 x 10“V 


+9.1969 X 10"V - 7.5711 x 10“V) exp(-a;) , 

' ZTT^ 

where x = E/T. 

Figure 2 presents a comparison of our electron-neutrino and neutrino-difference 


profiles with the estimate of ref. |jT^, which is 

E niE 


i^lye{E)/T = (6 X 10-") 


T V 4T 


^2 




(7.2) 


Our results exhibit two principle differences from those of ref. [^. First, the ana¬ 
lytical approximation ( |7.2| ) somewhat overestimates the electron-neutrino excess. In 
the canonical volume with a(t) = 2.0 MeV“^ at t = tf, we obtain a total electron- 
neutrino excess of about 2 x 10“^. In contrast, the analytical approximation (|7.2| ) 
gives about 4.8 x 10“^. Second, our distortions reach a maximum at a sightly lower 
energy. The peak in our neutrino-difference profile occurs at an energy value about 
25% below that of the peak in Eq. (fF^) . 

The total muon-neutrino excess we obtain in the canonical volume is about 1 x 
10“^ so that the total excess of electron neutrinos over muon neutrinos is also about 
1 X 10“^. Since the oscillation simulations in refs. |^, ^ were based on the 

initial production prohle ( |7.2|) , our current simulations lead to some modihcations of 
the earlier results. These effects are discussed in Sect IX. 

Figure 3 shows the time development of the nonequilibrium distortions. Various 
combinations of components of the average neutrino vector are plotted as functions 
of time. The curves begin at tg ^ 0.082 seconds and end at the hnal production 
time of tf ~ 1.312 seconds. The production proceeds smoothly and monotonically. 
The total number of excess electron neutrinos in the canonical volume corresponds 
to ((uo + "f^i)/2), while the total number of excess muon neutrinos orresponds to 
((uo ~ 'i^i) /2). The difference between electron and muon neutrinos is given by the 
component (vi), while their sum is (vq). The total excess of tau neutrinos is the 
same as for muon neutrinos, and hence is given by (^4) = ((uq — Ui) /2). Again, 
CP symmetry ensures that the results for antineutrinos are the same as those for 
neutrinos. 
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Overall, our results for pure production are similar to those of ref. [^. The 
hnal neutrino distortion productions in ref. were terminated at 0.1 MeV. When 
extrapolated to Tf, they are about 35% larger than ours but are otherwise comparable. 
As previously noted, the difference lies in the treatment of the electron mass. However, 
since we do not assume neutrinos are in thermal equilibrium at very early times, our 
early-time results differ. For instance, taking into account thermalizing neutrino 
interactions, we hnd that the T = 8 MeV electron-neutrino distortion is virtually 


zero, unlike the results in Figure 4 of ref. [|I^. For the same reason, we hnd a much 
smaller prohle at T = 4 MeV. Similarly, our results for Figures 5-8 of ref. produce 
curves dropping sharply for T > 4 MeV, without the 1/T^ behavior for early times. 


VIII. Results for Production with Oscillations 

This section discusses our numerical results for the production of nonequilibrium 
distributions when neutrino mixing is present. 

For A < 10“® eV^, our results resemble the 6 = 0 case discussed in Sect. VII. This 
follows because the vacuum term is small for tg < t < tf and A < 10“® eV^. The 
neutrino vectors, antineutrino vectors, and effective magnetic helds are all directed 
near the havor axis, i.e., the 1-axis. Since they point in a common direction, the 
cross-product terms v x By and w x By, in Eqs. ( b.8| ) and (|5.9|) are close to zero, and 
production proceeds as if there were no neutrino mixing. 

When A > 10“® eV^, some excess electron neutrinos convert to muon neutrinos 
during the production phase. As expected, the effect is greater for larger 6 . As a 
consequence of the neutrino conversion, the electron-neutrino and difference prohles 
are reduced while the muon-neutrino prohle is increased. 

Figure 4 displays the results for A = 10“® eV^ and sin^ 26 = 0.81. The curves 
representing the electron- and muon-neutrino sum prohle and the tan-neutrino prohle 
are almost identical to those in Figure 1. This also holds true at other values of A and 
6 . The reason is that the oscillation terms do not directly enter the expressions for 
these prohles, as can be seen from Eqs. (|5.8|) and (|5.9|) . Thus, oscillations ahect them 
only by feedback from the conversion among electron and muon neutrinos into the 
scattering processes in the Boltzmann equation. This secondary ehect is relatively 
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small. 

We have constructed an analytical approximation to the production profiles when 
0 < 6 ^ < 7 r/ 4 , given the pure-production results. It is based on the following idea. In a 
small time interval during production, a new excess of electron neutrinos is supplied 
to the system. This corresponds to adding small components to neutrino vectors 
along the 1-axis. These small additional components rotate around the associated 
magnetic helds. Since the production is continuous and since the time scale for 
oscillations is much shorter than the time scale for production, when averaged over 
an oscillation time the new small components lead to contributions aligned with the 
magnetic fields. Thus, ANME configurations are achieved. The averaging effect 
is equivalent to projecting the production vectors onto their magnetic fields. To 
construct analytical ANME configurations, we exploit the scheme presented in ref. 
T3| and mentioned in Sect. IV: the partial cancellation between the neutrino-neutrino 


and CP-asymmetric terms means that reasonably accurate ANME configurations can 
be obtained using only the vacuum and CP-symmetric terms. Then, for each E one 
projects the new production contributions onto the analytical ANME configurations. 
This must be done continuously during the production phase, and the results summed. 

It is useful to find a simpler approach based on the final pure-production profile. 
Then, one can employ the hts in Eq. ( [7.2[ ). In much of the parameter region of our 
simulations, the magnetic fields governing neutrino oscillations are changing. Initially, 
the magnetic fields point along the 1-axis. Hence, ANME configurations are flavor 
eigenstates. For A > 10“® V^, magnetic helds rotate toward A during the later 
stages of production. If most of the production occurs before this rotation, then the 
production prohle for 6 ^ 7 ^ 0 is obtained by rotating the pure-production case, for 
which 6 = 0, onto ANME conhgurations. If most of the production occurs after the 
rotation, then the production prohle for 6 ^ 0 is obtained by projecting the pure- 
production case onto ANME conhgurations. In general, early production is rotated 
onto ANME conhgurations, while later production is projected. Hence, we expect 
the true production prohle to lie between the projected and rotated curves. 

The analytical approximation to the data is implemented as follows. It is un¬ 
derstood that all vectors in this paragraph have three components. Dehne the CP- 
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symmetric part -B+ of the magnetic field by 


B+ {E,t) — ycp+{E)ei 


( 8 . 1 ) 


Let {E,tf) be the vector obtained in numerical simulations for pure production 
with 9 = 0. Then, the rotated vector fy {E,tf) is given by 




-B+ [E.tf) 


( 8 , 2 ) 


\E+ [E.tf) 

where the ± sign is the sign of B+ {E,tf). The projected vector Vp{E,tf) 

is given by 

77(0) . R, ^ 

(8.3) 


Vp (E, t,) = (b, t,) 


B+ ■ {E,tf) 

The analytical approximation based on rotation, with data approximated using 
Eq. (^.21) , works quite well for small A or small 6. For example, when A < 10“^ 
eV^, the agreement for vi {E, tf) and V 2 {E,tf) is at the one-percent level even for the 
large-angle case with sin^ 29 = 0.81. The same accuracy is obtained for the case with 
A = 10-6 eV^ with sin^ 29 = 0.25. 

Figures 5a and 5b show prohle plots of the components vi and V 2 for the case 
A = 10-6 eV^, sin^26* = 0.81. As expected, the data lie between the two analytical 
approximations based on rotation and on projection. In the region where F is a few 
MeV, the analytical rotation method overestimates the vectors by about 10%. The 
error reduces to about 5% in the energy range near 8 MeV, while for F > 11 MeV the 
method is accurate to at least one percent. The higher-energy vectors are reproduced 
better because their magnetic helds undergo relatively less rotation. 

The analytical approximation also works at earlier times during the production 
process. The idea is to use Eq. (|8.2|) or with the hnal time tf replaced by 

an arbitrary time t in the production interval tg < t < tf. For example, consider 
the time t ~ 0.5 seconds in the middle region of the production phase for the case 
with A = 10-"^ eV^ and sin^ 29 = 0.81. As expected, our data show that the values 
of Vi and V 2 again fall between the curves obtained from the rotation and projection 
approximations. For E > A MeV, the rotation method fits the data to about 10%. For 
F < 4 MeV, the rotation results are larger than the data, with signihcant deviations 
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when E < 2 MeV. For E <2 MeV, the projection method underestimates the data 
by about 30%. The projection approximation is better for low-energy neutrinos here 
because the corresponding ANME configurations undergo signihcant early rotation 
during production. 

In the parameter region with A > 10“^ eV^, the full simulation of the production 
phase involves prohibitive computer time. The production profile should he closer 
to the projection method based on Eq. (|8.3|). For A > 10“^ eV^, the projection 
approximation should be reasonably accurate. Indeed, in a test run for sin^ 29 = 0.81 
and A = 10“^ eV^, the data were reproduced within statistical uncertainties. 

The analytical approximation and the insights obtained in our earlier work, sum¬ 
marized in Sect. IV, suggest several predictions for neutrino behavior during the 
production phase. First, neutrinos should maintain themselves in ANME configura¬ 
tions. Second, the flavor development of neutrinos should be smooth. This follows 
from the ANME property and the slow evolution of the effective magnetic fields 
during the production time interval from tg to tf. Third, alignment should hold. 
This is a consequence of the ANME property and the approximate alignment of the 
effective magnetic helds. Fourth, there should be planarity. This follows because 
instantaneous nonlinear mass eigenstates have ^3 = 0 and so, since the neutrinos are 
in ANME conhgurations, ^3 should be much smaller than vi and V 2 - Finally, CP 
asymmetry should be suppressed. The ANME property implies the existence of the 
mechanism for CP-asymmetry cancellation given in ref. 0 . 

Data obtained during our production runs conhrm these hve predictions. All the 
data display smooth neutrino-flavor behavior. Figure 6 shows the components of 
the average neutrino vector as a function of time for the case with A = 10“® eV^ 
and sin^ 20 = 0.81. The smooth behavior is evident. The planarity feature is also 
evident since the third component is small. In fact, the third component is at least 
three orders of magnitude smaller than the first component. The antineutrino-vector 
components are omitted from Figure 6 because they are indistinguishable from the 
vector components on the scale of the figure. The difference between neutrino and 
antineutrino vectors is typically four or more orders of magnitude smaller than the 
vectors themselves. This reflects the CP-suppression mechanism. We also examined 
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data taken during the A = 10“^ eV^ and sin^ 26 = 0.81 production cycle, which 
is expected to exhibit the most sensitive dependence on the ANME and alignment 
properties. The data reveal that these properties hold well. The test results are 
similar to those shown in Table IV of ref. 0- 


IX. Results for Pure Oscillations after Production 

This section describes results for oscillations governed by Eq. ( p.2|) but based on 
the hnal production prohle obtained from the Boltzmann equation with oscillations. 
We also discuss some differences with results in our earlier works |^, [13 


m 


which the production prohle is treated in a step-function approximation based on the 
analytical approximation of ref. given by Eq. m- For purposes of comparison 
with the earlier work, we introduce in this section the normalized vectors r^{t) dehned 
by 


RP) = 


m) 


(9.1) 


imi) ■ 

In this equation, (|n(t)|) = Xj \ is the total number of neutrinos in the canonical 
comoving volume at time t, which is constant during the oscillation phase t > tf 
because neutrinos are neither created nor destroyed in the comoving volume. 

Figure 7 displays the components and r ^2 as a function of time for the case 
with A = 10“^^ eV^ and sin^ 29 = 0.81. The solid lines represent our new data, while 


the dashed lines represent data from earlier work |^. The hgure shows that a havor 
conversion occurs about 5 seconds earlier in our new data. 

The qualitative features of the two sets of data agree. In both sets of data, the 
behavior is evolutionary and smooth. At early times with t < 40 seconds, neutrinos 
are in approximate flavor eigenstates and hence the neutrino vectors he along the 
1-axis. The relatively small initial vacuum term grows with time and eventually 
dominates, so that for t > 120 seconds the neutrinos are instead in approximate 
vacuum-mass eigenstates with their vectors lying along A. Both sets of data display 
a feature called precocious rotation, i.e., the rotation to vacuum-mass eigenstates 
occurs signihcantly earlier than in the absence of the nonlinear term. The results 
differ signihcantly from data taken in the absence of neutrino-neutrino interactions. 
Without the nonlinear term, magnetic helds are dominated by the vacuum term only 
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for t > 500 seconds, damped oscillations occur, and neutrinos attain vacuum behavior 
only after t > 1500 seconds. 

The reason for the earlier flavor conversion in the new data is the use of our 
new production prohle. As noted in the discussion below Eq. (|7^) , there are two 
main differences between our new prohle and the one obtained in ref. [^: the total 
production is smaller, and the average energy is lower. See Figure 2. The change 
in the total production largely cancels in the ratio fy. It is the energy shift that 
is responsible for the 5-second time difference appearing in Figure 7. The partial 
cancellation between the terms lAj, and Vcp- means that the precocious rotation 
occurs when |A/2£'| and \Vcp+{E)\ become approximately equal When the 

average energy of the prohle is smaller, the vacuum term is larger and so the rotation 
to vacuum-mass eigenstates occurs sooner. 

A similar ehect appears for other values of A and sin^ 29. For example, the curves 
for Vyi and r ^2 in fhe case with A = 10“® eV^ and sin^ 29 = 0.81 are shifted relative 
to our earlier results by about one second in the crossover region between 3 and 16 
seconds. This is again due to the smaller average energy of the prohle, which causes an 
early rotation to vacuum-mass eigenstates because ANME conhgurations align sooner 
with A. Beyond 16 seconds, the two sets of data coincide. Likewise, the curves for r^i 
and r ^2 in the case with A = 10“® eV^ and sin^ 29 = 0.81 display time shifts of about 
0.1 seconds in the crossover region between 0.3 and 1.5 seconds. Beyond 2.0 seconds, 
results again coincide. For cases with smaller mixing angles, the time-diherence ehect 
is comparable but is less pronounced in plots because the neutrino havor changes less. 
We remark in passing that the ehect is absent for the case with A = 10“^ eV^ and 
sin^ 29 = 0.81 because the overlap between the new data and those of ref. occurs 


when neutrinos are already near vacuum-mass eigenstates. 


The results we have obtained conhrm all the qualitative conclusions of ref. [13 


summarized in Sect. IV. The numerical results of m are accurate, except for the 
above-mentioned time-shift ehect in the crossover region. 


The analytical approximations obtained in the present work and in ref. |13| pro¬ 


vide a means to determine neutrino vectors to an accuracy better than the intrinsic 
methodological uncertainties discussed in Sect. II. There are two stages involved. 























First, given a value of A and sin^26', the neutrino-distortion production profile can 
be determined via the techniques provided in Sect. VIII. Second, results for post¬ 
production times t > tf can be found using the method of Sect. VI.E in ref. ||13|| . 
This method gives the neutrino vectors as 

. ( 9 . 2 ) 

|A+ IE, tj| 


where 5+ {E,t) is dehned in Eq. (p.l|) , and where the overall sign is that of the dot 
product v{E,tf) ■ {E,tf). The antineutrino vectors are set equal to the neutrino 

vectors, w{E,t) = v{E,t), thus guaranteeing perfect CP-suppression. The third 
components of vectors are zero, thus ensuring planarity. The second stage of the 
method reproduces numerical results to about two decimal places []^. Via this two- 
step procedure, reasonably accurate numbers for all quantities can be obtained. 


X. Summary 

In this work, we have treated the nonequilibrium thermodynamic statistics of neu¬ 
trinos in the early Universe in the presence of neutrino mixing and oscillations. Differ¬ 
ential equations were derived that determine the deviations from standard equilibrium 
statistics. They are Eqs. (|5.8|) and (|5.9| ). We have resorted to numerical methods to 
solve these equations, since they are complicated and nonlinear. We summarize here 
the behavior of the solutions and our methods for approximating them analytically. 

Soon after the start of the production of nonequilibrium distortions, neutrinos and 
antineutrinos achieve configurations that are approximate nonlinear mass eigenstates 
(ANME). As the Universe expands, various interaction terms change. However, the 
neutrinos remain in ANME configurations during these changes, so the flavor content 
evolves smoothly. 

For A < 10“® eV^, the production of nonequilibrium distortions occurs when the 
vacuum term is relatively small. This implies that the production phase proceeds as 
in the case without neutrino mixing, producing neutrinos close to flavor eignenstates. 
During the oscillation phase, the vacuum term increases while the interaction effects 
decrease. Eventually, the vacuum term dominates. Neutrinos therefore gradually 
evolve from flavor eigenstates to vacuum-mass eigenstates. For A < 10“® eV^, the 
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rotation occurs earlier than one might expect if neutrino-neutrino interactions are 


neglected. This precocious rotation, observed and explained in ref. |]^, is conhrmed 
in our simulations. 

For A > 10“® eV^ oscillations play a signihcant role already during the production 
phase. Since the nonequilibrium distortions are induced by the weak interactions, 
the excess neutrinos are generated as flavor eigenstates. However, the time scale 
for oscillations is much shorter than that for production, so the neutrinos quickly 
oscillate, decohere, and achieve ANME conhgurations. For A < 10“^ eV^ much of the 
production occurs while ANME conhgurations point in the havor direction. Hence, a 
good approximation to the hnal production prohle in this case is to rotate the hnal 
pure-production (zero mixing) results onto ANME conhgurations. For A > 10“^ 
eV^ most of the production occurs while ANME conhgurations point in the mass- 
eigenstate direction. Hence, a good approximation to the hnal production prohle 
in this case is to project the pure-production results onto ANME conhgurations. In 
the region 10“^ eV^ < A < 10“^ eV^ these two approximations straddle the true 
production results. 

The above approximations to the hnal production prohle can be made analytical 
by using the cancellation mechanism of ref. 0 so that mass eigenstates are con¬ 
structed based solely on the vacuum term and the CP-conserving interaction between 
the charged leptons and the neutrinos. For the parameter region of our simulations, 
we hnd that the rotation method reproduces the data extremely well for A < 10“^ 
eV^. For A 10“® eV^ and sir? 26 = 0.81, it overestimates results by at most 10%. 
One can compensate for this ehect. We have also described in Section VHI how to 
extend the analytical method to the production-oscillation results for larger A. 

During the oscillation phase for A > 10“® eV^, neutrinos continue to rotate to¬ 
ward vacuum-mass eigenstates. The behavior therefore continues to be smooth and 
evolutionary. For A > 10“^ eV^, we expect the production phase to render neutri¬ 
nos close to vacuum-mass eigenstates. Since the vacuum term dominates thereafter, 
neutrinos remain in vacuum-mass eigenstates during the oscillation phase. The result 
should be constant flavor behavior. Indeed, short test simulations conhrm this. 

The improved treatment of the production phase means that our results for un- 
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normalized quantities are now accurate to within systematic uncertainties. This rep¬ 
resents an improvement in absolute accuracy of a factor of about twenty. 

The approximate analytical method of ref. can be used to obtain to good pre¬ 
cision the results for the oscillation phase. By combining this method with the new 
analytical treatment of the production, we obtain a complete analytical approxima¬ 
tion for the neutrino and antineutrino vectors over the entire time period of interest. 

All the qualitative neutrino flavor behavior of refs. 


has been conhrmed. 


Features including smooth evolutionary behavior, alignment, planarity, maintenance 
of ANME configurations, and CP suppression are observed not only in the oscillation 
phase but also during the production phase. 


We have uncovered one quantitative difference with ref. neutrinos rotate to 
vacuum-mass eigenstates about 5% earlier in time. The reason is that the approximate 


initial distortion distributions used in ref. are shifted somewhat to higher energies 
than those produced in our present numerical data. This caused the vacuum term to 
dominate somewhat later in the simulations of ref. []^. Asymptotically, the results 
of our current and previous work agree. 

Within the parameter region of this paper, we have found little generation of CP 
asymmetry. This means that nonequilibrium distortions are unlikely to affect nucle¬ 
osynthesis substantially. Although one might a priori expect neutrino oscillations to 
have observable effects, the CP-suppression mechanism we uncovered in our simula¬ 
tions implies that the impact on the primordial helium abundance is likely to be at 
the same level as found in ref. [|T^, where neutrino mixing was absent. 
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Appendix A: Coefficients for Three-Flavor Pure Production 

This appendix provides the A, B, and C coefficients appearing in Eq. They 

are fnnctions of the scaled initial and final nentrino energies ^ and This fnnctional 
dependence is dropped here to save space. In what follows, it is nsefnl to define the 
electroweak conpling constants 

a = (1 + sin^ 6'u,)^ , 6 = (2 sin^ 0^,)^ , c = (1 — sin^ 6*^)^ , (A.l) 

where the sqnared weak-mixing angle is 

sin^ ^ 0.2325 . (A.2) 


The A coefficients are given by 


^ ^ p(a + fc) + 17 \g 

*^e \ 3 / T 


^ - ( Hb + c) +17 \ E 


(A.3) 


The B coefficients are 


E 


, E fllE \ 

Bp, = B^, = {a + b) — - Ij exp ^ ^ 


E /IIE \ / E 

Bp^ = B^^ = Bp^ = B^^ = {b + c) — - - ij exp 


(A.4) 


The C coefficient are expressed in terms of functions Qi, presented in Eq. ( |A.12| ) 
below. For the electron neutrino, we find 

= CIqi + Cg(73 + 2Cg(72 + 3^0 , 

Ci/ei/,, =91 + 2(73 + 2^2 + 3^0 5 Cj^ei^T = 9i + 2^3 + 2(72 + 3(7o , 

Cu,p, = Agi + CIqq , Cy,p^ = 2^1 + (7o , Cy,p^ = 2^i + 770 . (A.5) 

For the muon neutrino, we find 


“ ^ 1.31 + C^fi'3 + 2C^fi'2 + 35^0 , 
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— S'! + 2(73 + 2(72 + 3(7o , — gi + 2773 + 2(72 + 3(7o , 

= 2^1 + (7o , = 4^1 + 1 C'l/^.Pr = 2^1 + % ■ (A.6) 

For the tau neutrino, we find 

= C^gi + C^gs + 2C^^2 + ^go , 

= gi + 2^3 + 2(72 + 3(7 o , ~ + 2(73 + 2^2 + 3^o > 

Cl^rPe = 2(71 + % , C'l^rP^ = 2(71 + ^0 5 C'^rPr = 4^1 + C'^^O • (A-T) 

For the electron antineutrino, we find 

Cpel^e = 4^1 + Clgo , C'Pel'M ^ ^ ^0 , Cpet^T = 2(7l + ^0 , 

CpePe = + (Fg (73 + 2Cg (72 + 3^0 5 

C'pepM = + 2^3 + 2^2 + 3(7 o , C'pePr = + 2^3 + 2^2 + 3(7o . (A.8) 

For the muon antineutrino, we find 

Cp^i^e = 2(7 i + ^0 , = 4^1 + C^^go , Cy^y^ = 2gi + % , 

Cpppp = + ^*^^3 + 2C^(72 + 3^0 , 

Cp^Pe = + 2^3 + 2^2 + 3(7 o , C'p^p^ = gi + 2^3 + 2(72 + ?>% . (A.9) 

For the tau antineutrino, we find 


Cp^i/e — 2(7 i + (7o 5 Cy^y^ — 2(7i + (7o , Cy^y^ — Aq I + C^(7o , 

Cy^y^ = C^gi + C^gs + 2C^^2 + ^go , 

Cy^y^ = gi + 2^3 + 2^2 + 3^0 , Cy^y^ = + 2^3 + 2^2 + 3^0 • (A. 10) 

In Eqs. ( |A.5D -( |A.10| ), we use the abbreviations 

C*g = ( 2 + 6+6 , C*^ = 6+C+6 , C*g = 2(a+6) + 10 , = 2(6+c) + 10 . (A.11) 


The g are scaled versions of the functions g of ref. |15]: 


go 


m)-rsm)M 4 
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(A.12) 



For completeness we give the unbarred 


gi (f, f)=f ~ y‘^)‘ 

^ ^ j\v\ ^ ' 


g2 


(f. f) = - /” (2“ + + »') . 


^3 


(hf) - 


l\v\ Ay 


X + [w"^ + Aw — ^ + 2 (w — 1) y'^ + y'^^ (A-13) 


where the abbreviations 




(-B + E’) 

w = 


E 

X = — 

T 


(A. 14) 


rjn 7 rji 7 

are used. The lower limit of each integration in Eq. ( [A.13|) is the absolute value of v. 
Note that the above expression for g'^ differs from that in ref. |M, which is incorrectly 
typeset. 


Appendix B: Coefficients for Combined Production-Oscillation 

Here, we express the coefficients Aab, A^i, Ba, Ba, Cab, Cai, Cab, and Cai in the 
vector basis in terms of the coefficients and Cy^v^,. Since all coefficients 

are functions of the scaled initial and final neutrino energies ^ and we drop this 
dependence everywhere to save space. 

The nonzero Aab and are 



(BA) 


34 













and 


1 1 

^00 = ’ ^oi ^ ^io = ^ (^Pe ~ ^Pm) ’ ^23 = 


The nonzero Br, and Bn are 


(B.2) 


Bq — By^ + Bjj 


•fx ) 


Bl — By^ — By 


Bi — By^ 


(B.3) 


and 


Bq — By^ + By^ , Bl — By^ “ By^ , B^ — B y^ 
The remaining nonzero coefficients for nentrinos are 

1 


Coo — 


^01 = 


{Cy^y^ + Cy^y^ T Cy^y^ T C y ^y ^ , 

(Cy^y^ + Cy^y^ — Cy^y^ — C y ^y ^ , 

ClO = - {Cy^y^ — Cy^y^ + Cy^y^ — C y ^y ^ , 

Cll = - {Cy^y^ — Cy^y^ — Cy^y^ T Cy^y^ , 

C04 = + Cy^y^^ , Cl4 = (Cy^y^ — C y ^y^ , 


(B.4) 


C4O — 2 (Cj/eJ^T + Cy^y^ , C4I — - [Cy^y^ — Cy^y^ , C44 — (T, 

CqO 2 (Cj^ePe T Cj/^pg + Cy^y^ + J 


V-tVt 1 


^1 = 


(Cy^y^ + Cy^y^ Cy^y^ Cy^y^ , 

ClO = 2 {Cy^y^ — Cy^y^ + Cl/gP^ — C y ^y ^ , 


C'n = 


Cy^y^ Cy^y^ + Cy^y^ , 

C04 = (py^y^ + Cy^P , Cl4 = (Cl/gP^ — Cy^yj^ , 


C40 = ^ (agp. + a,p.) , C41 = ^ (agp. - a,p.) , C 44 = a 


Vr^T 


(B.5) 
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For antineutrinos, the remaining nonzero coefficients are 


C^oo = 


{Cu^Us + + Cy^y^ + C y ^y ^ 


^01 = 2 — Cy^y^ — Cy^y^ 

C*l0 = 2 ~ ~ 

Cll = - — Cy^y^ — Cy^y^ + C y ^y ^ 


Cu = (c,^u.+C,^u?j , ,Cu={q 


— r'- 

VeVr 


C4O — - (Cy^y^ + Cy^y^^ , , C4I — “ (Cy^y^ — Cy^y^^j , , C44 — Cy 


Coo — 2 (CpePe + Cp^Pe + Cp^Pp + Cp^Pp) 


Coi = 


C'io = 


^11 = 


(c ntiePe ffi Cp^Pg CpgP^ Cp^p^^ 
^ (CpePe ~ Cp^Pg + CpgP^ — Cy^y^ 

^ (CpePe ~ Cp^Pg — Cy^y^ + Cy^y^ 


C 4 O — 


Co4 — (Cp^p^ + Cy^y^^ , , Cl4 — (Cy^y^ — Cp^p^^ , 

i (&.P, + Cp,p,) , .Cji = i(C,.,,-4p,) , ,C4,=0 


(B.6) 
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FIGURE CAPTIONS 


Figure 1: Pure-production profiles. For zero mixing angle and at T = 0.75 MeV, the 
excess-neutrino number per unit energy in the canonical comoving volume is shown 
as a function of energy. The solid line is electron neutrinos, the dashed line is muon 
neutrinos, the dotted line is their difference, and the dashed-dotted line is their sum. 


Figure 2: Comparison of pure-production prohles. For zero mixing angle and at 
T = 0.75 MeV, the excess-neutrino number per unit energy in the canonical comoving 
volume is shown as a function of energy. The solid line is our electron-neutrino results, 
the dotted line is the difference between our electron- and muon-neutrino results, and 


the dashed line is the analytical approximation of ref. [|l6 


Figure 3: Time development of nonequilibrium distortions. The solid line is the total 
excess of electron neutrinos in the canonical comoving volume as a function of time. 
The dashed line is the total muon-neutrino excess, the dotted line is the difference 
between electron and muon neutrinos, and the dashed-dotted line is their sum. 

Figure 4: Production prohles for A = 10“® eV^ and sin^ 29 = 0.81. The excess- 
neutrino number per unit energy in the canonical comoving volume is shown as a 
function of energy. The solid line is electron neutrinos, the dashed line is muon 
neutrinos, the dotted line is their difference, the dashed-dotted line is their sum, and 
the short-dashed line is tau neutrinos. 


Figure 5: Analytical approximations to prohle plots for A = 10 ® eV^, sin^ 29 = 0.81. 
The solid lines are the data, the dashed lines are the projection approximation 


and the dotted lines are the rotation approximation (|8.2|) . (a) The hrst component 
of the neutrino vector as a function of energy, (b) The second component of the 
neutrino vector as a function of energy. 


Figure 6: Components of the average neutrino vector as a function of time for the 
case with A = 10“® eV^ and sm^29 = 0.81. The dashed-dotted line is the zeroth 


38 






component, the dotted line is the hrst component, the solid line is the second com¬ 
ponent, and the short-dashed line is the fourth component. The third component is 
indistinguishable from zero on the scale of the hgure. 


Figure 7: Components of the normalized average neutrino vector f^{t) as a function of 
time for the case with A = 10“^^ eV^ and sin^ 26 = 0.81. The solid lines represent the 


new data, while the dashed lines represent data from earlier work |^. The curves 
above the horizonal axis are the component r^i, while the curves below it are the 
component r^ 2 - 
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